#include <cmath>
#ifdef _MPI
#include <mpi.h>
#endif
#include "cgnslib.h"

#include "predefined.hpp"
#include "global.hpp"
#include "node.hpp"
#include "mesh.hpp"
#include "preprocess.hpp"
#include "solution.hpp"
#include "flux.hpp"



using namespace std;

Preprocess::Preprocess(){
    m_scheme(s);
    m_msg.reserve(MAXCHAR);
}

Preprocess::~Preprocess(){
    m_msg="Preprocess destructor function has been called.";
    logCtrl.info(m_msg.c_str());
}



void Preprocess::m_scheme(Scheme &s){
    if (STENCIL>2){
    s.d[0][2][0] = -3.0/2.0;
    s.d[0][2][1] =  2.0;
    s.d[0][2][2] = -1.0/2.0;
    s.d[1][1][0] = -1.0/2.0;
    s.d[1][1][1] =  0.0;
    s.d[1][1][2] = -s.d[1][1][0];
    s.d[2][0][0] = -s.d[0][2][2];
    s.d[2][0][1] = -s.d[0][2][1];
    s.d[2][0][2] = -s.d[0][2][0];
    }

    //4 points
    if (STENCIL>3){
    s.d[0][3][0] = - 11.0/6.0;
    s.d[0][3][1] =   18.0/6.0;
    s.d[0][3][2] = -  9.0/6.0;
    s.d[0][3][3] =    2.0/6.0;
    s.d[3][0][0] = - s.d[0][3][3];
    s.d[3][0][1] = - s.d[0][3][2];
    s.d[3][0][2] = - s.d[0][3][1];
    s.d[3][0][3] = - s.d[0][3][0];

    s.d[1][2][0] = - 2.0/6.0;
    s.d[1][2][1] = - 3.0/6.0;
    s.d[1][2][2] =   6.0/6.0;
    s.d[1][2][3] = - 1.0/6.0;
    s.d[2][1][0] = - s.d[1][2][3];
    s.d[2][1][1] = - s.d[1][2][2];
    s.d[2][1][2] = - s.d[1][2][1];
    s.d[2][1][3] = - s.d[1][2][0];
    }

    //5 points
    if (STENCIL>4){
    s.d[0][4][0] = -25.0/12.0;
    s.d[0][4][1] =  4.0;
    s.d[0][4][2] = -3.0;
    s.d[0][4][3] =  16.0/12.0;
    s.d[0][4][4] = -3.0/12.0;
    s.d[1][3][0] = -3.0/12.0;
    s.d[1][3][1] = -10.0/12.0;
    s.d[1][3][2] =  18.0/12.0;
    s.d[1][3][3] = -0.5;
    s.d[1][3][4] =  1.0/12.0;
    s.d[2][2][0] =  1.0/12.0;
    s.d[2][2][1] = -8.0/12.0;
    s.d[2][2][2] =  0.0;
    s.d[2][2][3] = -s.d[2][2][1];
    s.d[2][2][4] = -s.d[2][2][0];
    s.d[3][1][0] = -s.d[1][3][4];
    s.d[3][1][1] = -s.d[1][3][3];
    s.d[3][1][2] = -s.d[1][3][2];
    s.d[3][1][3] = -s.d[1][3][1];
    s.d[3][1][4] = -s.d[1][3][0];
    s.d[4][0][0] = -s.d[0][4][4];
    s.d[4][0][1] = -s.d[0][4][3];
    s.d[4][0][2] = -s.d[0][4][2];
    s.d[4][0][3] = -s.d[0][4][1];
    s.d[4][0][4] = -s.d[0][4][0];
    }

    //7 points
    if (STENCIL>6){
    s.d[0][6][0] = - 2.192280339;
    s.d[0][6][1] =   4.748611401;
    s.d[0][6][2] = - 5.108851915;
    s.d[0][6][3] =   4.461567104;
    s.d[0][6][4] = - 2.833498741;
    s.d[0][6][5] =   1.128328861;
    s.d[0][6][6] = - 0.203876371;
    s.d[1][5][0] = - 0.209337622;
    s.d[1][5][1] = - 1.084875676;
    s.d[1][5][2] =   2.147776050;
    s.d[1][5][3] = - 1.388928322;
    s.d[1][5][4] =   0.768949766;
    s.d[1][5][5] = - 0.281814650;
    s.d[1][5][6] =   0.048230454;
    s.d[2][4][0] =   0.049041958;
    s.d[2][4][1] = - 0.468840357;
    s.d[2][4][2] = - 0.474760914;
    s.d[2][4][3] =   1.273274737;
    s.d[2][4][4] = - 0.518484526;
    s.d[2][4][5] =   0.166138533;
    s.d[2][4][6] = - 0.026369431;
    s.d[3][3][0] = - 0.0208431427703117643;
    s.d[3][3][1] =   0.166705904414580469;
    s.d[3][3][2] = - 0.770882380518225552;
    s.d[3][3][3] =   0.0;
    #ifdef _OPTDRP
    s.d[0][6][0] = - 2.2219902954809613859765704546589404;
    s.d[0][6][1] =   4.7937139708231946144678659038618207;
    s.d[0][6][2] = - 4.8887154219014803402387769892811775;
    s.d[0][6][3] =   3.7241945556599418409859936218708754;
    s.d[0][6][4] = - 1.9475764115885354499368986580520868;
    s.d[0][6][5] =   0.64080276257284296548277779947966337;
    s.d[0][6][6] = - 0.10042916008504576552695652935653925;
    s.d[1][5][0] = - 0.19275393397114015670013031922280788;
    s.d[1][5][1] = - 1.1617462897578052505309642583597451;
    s.d[1][5][2] =   2.2833737916894651931443149806000292;
    s.d[1][5][3] = - 1.4942869230903372823604513541795313;
    s.d[1][5][4] =   0.79138992627937021495654335012659431;
    s.d[1][5][5] = - 0.2681591974297292679807469539809972;
    s.d[1][5][6] =   0.042182626280174218003082842187723145;
    s.d[2][4][0] =   0.043961792318923609279579522990388796;
    s.d[2][4][1] = - 0.44665643971666746381288248812779784;
    s.d[2][4][2] = - 0.50947801953384086459664104040712118;
    s.d[2][4][3] =   1.2919072955902599630917393369600177;
    s.d[2][4][4] = - 0.51171625718486879463853256311267614;
    s.d[2][4][5] =   0.1551341504041561236704183102119714;
    s.d[2][4][6] = - 0.02315252187794980542889788921456784;
    s.d[3][3][0] = - 0.019408956520410941487537255049161731932458;
    s.d[3][3][1] =   0.16096915941497709928348235352998026106316;
    s.d[3][3][2] = - 0.76371144926872137410435294191247532632895;
    s.d[3][3][3] =   0.0;
    #endif
    s.d[3][3][4] = - s.d[3][3][2];
    s.d[3][3][5] = - s.d[3][3][1];
    s.d[3][3][6] = - s.d[3][3][0];
    s.d[4][2][0] = - s.d[2][4][6];
    s.d[4][2][1] = - s.d[2][4][5];
    s.d[4][2][2] = - s.d[2][4][4];
    s.d[4][2][3] = - s.d[2][4][3];
    s.d[4][2][4] = - s.d[2][4][2];
    s.d[4][2][5] = - s.d[2][4][1];
    s.d[4][2][6] = - s.d[2][4][0];
    s.d[5][1][0] = - s.d[1][5][6];
    s.d[5][1][1] = - s.d[1][5][5];
    s.d[5][1][2] = - s.d[1][5][4];
    s.d[5][1][3] = - s.d[1][5][3];
    s.d[5][1][4] = - s.d[1][5][2];
    s.d[5][1][5] = - s.d[1][5][1];
    s.d[5][1][6] = - s.d[1][5][0];
    s.d[6][0][0] = - s.d[0][6][6];
    s.d[6][0][1] = - s.d[0][6][5];
    s.d[6][0][2] = - s.d[0][6][4];
    s.d[6][0][3] = - s.d[0][6][3];
    s.d[6][0][4] = - s.d[0][6][2];
    s.d[6][0][5] = - s.d[0][6][1];
    s.d[6][0][6] = - s.d[0][6][0];
    

    }

    #ifdef _RK46
    //RK46
    //4 steps
    s.rk[0][0][0] = 0.0;
    s.rk[0][0][1] = 0.5;
    s.rk[0][0][2] = 0.5;
    s.rk[0][0][3] = 1.0;
    s.rk[0][1][0] = 1.0 / 6.0;
    s.rk[0][1][1] = 1.0 / 3.0;
    s.rk[0][1][2] = 1.0 / 3.0;
    s.rk[0][1][3] = 1.0 / 6.0;
    s.rk[0][2][0] = 0;
    s.rk[0][2][1] = 0;
    s.rk[0][2][2] = 0;
    s.rk[0][2][3] = 0;
    //6 steps
    s.rk[1][0][0] = 0.0;
    s.rk[1][0][1] = 0.353323;
    s.rk[1][0][2] = 0.999597;
    s.rk[1][0][3] = 0.152188;
    s.rk[1][0][4] = 0.534216;
    s.rk[1][0][5] = 0.603907;
    s.rk[1][1][0] = 0.0467621;
    s.rk[1][1][1] = 0.137286;
    s.rk[1][1][2] = 0.170975;
    s.rk[1][1][3] = 0.197572;
    s.rk[1][1][4] = 0.282263;
    s.rk[1][1][5] = 0.165142;
    s.rk[1][2][0] = 0;
    s.rk[1][2][1] = 0;
    s.rk[1][2][2] = 0;
    s.rk[1][2][3] = 0;
    s.rk[1][2][4] = 0;
    s.rk[1][2][5] = 0;
    #else
    //rk56
    //5 steps
    s.rk[0][0][0] = 0.0;
    s.rk[0][0][1] = -0.6051226;
    s.rk[0][0][2] = -2.0437564;
    s.rk[0][0][3] = -0.7406999;
    s.rk[0][0][4] = -4.4231765;
    s.rk[0][1][0] =  0.2687454;
    s.rk[0][1][1] =  0.8014706;
    s.rk[0][1][2] =  0.5051570;
    s.rk[0][1][3] =  0.5623568;
    s.rk[0][1][4] =  0.0590065;
    s.rk[0][2][0] = 0.0;
    s.rk[0][2][1] = 0.2687454;
    s.rk[0][2][2] = 0.5852280;
    s.rk[0][2][3] = 0.6827066;
    s.rk[0][2][4] = 1.1646854;
    //6 steps
    s.rk[1][0][0] = 0.0;
    s.rk[1][0][1] = -0.4412737;
    s.rk[1][0][2] = -1.0739820;
    s.rk[1][0][3] = -1.7063570;
    s.rk[1][0][4] = -2.7979293;
    s.rk[1][0][5] = -4.0913537;
    s.rk[1][1][0] = 0.1158488;
    s.rk[1][1][1] = 0.3728769;
    s.rk[1][1][2] = 0.7379536;
    s.rk[1][1][3] = 0.5798110;
    s.rk[1][1][4] = 1.0312849;
    s.rk[1][1][5] = 0.15;
    s.rk[1][2][0] = 0.0;
    s.rk[1][2][1] = 0.1158485;
    s.rk[1][2][2] = 0.3241850;
    s.rk[1][2][3] = 0.6193208;
    s.rk[1][2][4] = 0.8034472;
    s.rk[1][2][5] = 0.9184166;
    #endif
    s.f[0][0][0] =  0.0;
    s.f[1][1][0] = -0.25;
    s.f[1][1][1] =  0.5;
    s.f[1][1][2] =  s.f[1][1][0];
    s.f[2][2][0] =  0.0625;
    s.f[2][2][1] = -0.25;
    s.f[2][2][2] =  0.375;
    s.f[2][2][3] =  s.f[2][2][1];
    s.f[2][2][4] =  s.f[2][2][0];
    if (STENFIL>5){
    s.f[3][3][0] = - 0.0238530481912;
    s.f[3][3][1] =   0.1063035787698;
    s.f[3][3][2] = - 0.2261469518087;
    s.f[3][3][3] =   0.2873928424602;
    s.f[3][3][4] =   s.f[3][3][2];
    s.f[3][3][5] =   s.f[3][3][1];
    s.f[3][3][6] =   s.f[3][3][0];
    }
    if (STENFIL>7){
    s.f[4][4][0] =   0.008228661760;
    s.f[4][4][1] = - 0.045211119360;
    s.f[4][4][2] =   0.120007591680;
    s.f[4][4][3] = - 0.204788880640;
    s.f[4][4][4] =   0.243527493120;
    s.f[4][4][5] =   s.f[4][4][3];
    s.f[4][4][6] =   s.f[4][4][2];
    s.f[4][4][7] =   s.f[4][4][1];
    s.f[4][4][8] =   s.f[4][4][0];
    }
    if (STENFIL>9){
    s.f[5][5][0] = - 0.002999540835;
    s.f[5][5][1] =   0.018721609157;
    s.f[5][5][2] = - 0.059227575576;
    s.f[5][5][3] =   0.123755948787;
    s.f[5][5][4] = - 0.187772883589;
    s.f[5][5][5] =   0.215044884112;
    s.f[5][5][6] =   s.f[5][5][4];
    s.f[5][5][7] =   s.f[5][5][3];
    s.f[5][5][8] =   s.f[5][5][2];
    s.f[5][5][9] =   s.f[5][5][1];
    s.f[5][5][10]=   s.f[5][5][0];
    }
}

void Preprocess::initialize(Mesh &msh, Node *node,Solution *sol, Flux *flux){
    int iLNode;
    int iDim,iEq,iVar;
    double ke,dis;

    //Initialize Q, lastQ, FRK, WRK with 0.
    for (iLNode=0;iLNode<msh.nLNode;iLNode++){
        for (iVar=0;iVar<NVAR;iVar++){
            (sol+iLNode)->Q[iVar]=0;
            (sol+iLNode)->lastQ[iVar]=0;
        }
        for (iEq=0;iEq<NEQ;iEq++){
            (flux+iLNode)->FRK[iEq]=0;
            (flux+iLNode)->WRK[iEq]=0;
        }
    }

    //Initialize the source
    
    for (iLNode=0;iLNode<msh.nLNode;iLNode++){
        if (inputCtrl.srcTyp==0){
            dis=0;
            for (iDim=0;iDim<NDIM;iDim++) dis+=pow(abs((node+iLNode)->xyz[iDim]-inputCtrl.srcPointPos[iDim]),2);
            dis=sqrt(dis);
            (sol+iLNode)->Q[DENSITY]+=inputCtrl.srcPamp*exp(-log(2.0)*(dis/inputCtrl.srcRadius)*(dis/inputCtrl.srcRadius));
            (sol+iLNode)->Q[PRESSURE]+=inputCtrl.srcPamp/inputCtrl.gamma*exp(-log(2.0)*(dis/inputCtrl.srcRadius)*(dis/inputCtrl.srcRadius));
        }
    }
}

void Preprocess::calJacobi(Mesh &msh, Node *node){
    int iDim,iZone,iLZone,iZNode,iLNode,nNode,dIdxNode;
    int iSten,iType,size[NDIM],i,j;
    int ijk[NDIM];
    double dl[NDIM],dlMax=TINYNUM,dlGlMax=TINYNUM,dlMin=HUGENUM,dlGlMin=HUGENUM;
    double invJ[NDIM][NDIM],Mat[DIMSQUR],J;
    //invJ | dx/dXi dx/dEta dx/dZeta |
    //     | dy/dXi dy/dEta dy/dZeta |
    //     | dz/dXi dz/dEta dz/dZeta |
    //invJ | dx/dXi dx/dEta |
    //     | dy/dXi dy/dEta |
    for (iLZone=0;iLZone<msh.nLZone;++iLZone){
        iZone=msh.lZoneList[iLZone];
        nNode=msh.nNode[iZone];
        //printf("Point B rank %d iZone %d\n",msh.cId,iZone);
        for (iDim=0;iDim<NDIM;++iDim) size[iDim]=msh.size[iZone][iDim];
        for (iZNode=0;iZNode<nNode;++iZNode){
            iLNode=msh.idxLNode[iZone]+iZNode;
            msh.calijk(iZone,iZNode,ijk);
            for (i=0;i<NDIM;++i){
                for (j=0;j<NDIM;++j){
                    invJ[i][j]=0;
                }
            }
            dIdxNode=1;
            for (iDim=0;iDim<NDIM;++iDim){
                iType=(node+iLNode)->locTyp[iDim];
                for (iSten=0;iSten<STENCIL;iSten++){
                    invJ[0][iDim]+=s.d[iType][STENCIL-iType-1][iSten]*(node+iLNode-(iType-iSten)*dIdxNode)->xyz[0];
                    invJ[1][iDim]+=s.d[iType][STENCIL-iType-1][iSten]*(node+iLNode-(iType-iSten)*dIdxNode)->xyz[1];
                    #ifndef _2D
                    invJ[2][iDim]+=s.d[iType][STENCIL-iType-1][iSten]*(node+iLNode-(iType-iSten)*dIdxNode)->xyz[2];
                    #endif
                }
                dIdxNode*=size[iDim];
            }
            //Metric
            for (i=0;i<NDIM;++i){
               for (j=0;j<NDIM;++j){
                   Mat[i*NDIM+j]=invJ[i][j];
               }
            }
            //if(iZone==0)printf("Point D: iZone %d J %lf z %lf\n",iZone,Mat[0],(node+iZNode)->xy[2]);
            J=math.determinant(Mat,NDIM);
            if (J==0){
                J=TINYNUM;
#               ifndef _2D
                m_msg="Zero Jacobian at: ("+to_string(ijk[0])+","+to_string(ijk[1])+","+to_string(ijk[2]);
                logCtrl.warning(m_msg.c_str());
#               else
                m_msg="Zero Jacobian at: ("+to_string(ijk[0])+","+to_string(ijk[1]);
                logCtrl.warning(m_msg.c_str());
#               endif
            }
            (node+iLNode)->J=J;
            //if(iZone==0)printf("Point E: iZone %d J %lf %d %d %d\n",iZone,(node+iZNode)->J,ijk[0],ijk[1],ijk[2]);
            //jac  |   dXi/dx   dXi/dy   dXi/dz |
            //     |  dEta/dx  dEta/dy  dEta/dz |
            //     | dZeta/dx dZeta/dy dZeta/dz |
            //jac  |   dXi/dx   dXi/dy |
            //     |  dEta/dx  dEta/dy |
            #ifndef _2D
            (node+iLNode)->ja[0][0]=( invJ[1][1]*invJ[2][2]-invJ[1][2]*invJ[2][1])/J;
            (node+iLNode)->ja[0][1]=(-invJ[0][1]*invJ[2][2]+invJ[0][2]*invJ[2][1])/J;
            (node+iLNode)->ja[0][2]=( invJ[0][1]*invJ[1][2]-invJ[0][2]*invJ[1][1])/J;
            (node+iLNode)->ja[1][0]=(-invJ[1][0]*invJ[2][2]+invJ[1][2]*invJ[2][0])/J;
            (node+iLNode)->ja[1][1]=( invJ[0][0]*invJ[2][2]-invJ[0][2]*invJ[2][0])/J;
            (node+iLNode)->ja[1][2]=(-invJ[0][0]*invJ[1][2]-invJ[0][2]*invJ[1][0])/J;
            (node+iLNode)->ja[2][0]=( invJ[1][0]*invJ[2][1]-invJ[1][1]*invJ[2][0])/J;
            (node+iLNode)->ja[2][1]=(-invJ[0][0]*invJ[2][1]+invJ[0][1]*invJ[2][0])/J;
            (node+iLNode)->ja[2][2]=( invJ[0][0]*invJ[1][1]-invJ[0][1]*invJ[1][0])/J;
            #else
            (node+iLNode)->ja[0][0]=( invJ[1][1])/J;
            (node+iLNode)->ja[0][1]=(-invJ[0][1])/J;
            (node+iLNode)->ja[1][0]=(-invJ[1][0])/J;
            (node+iLNode)->ja[1][1]=( invJ[0][0])/J;
            #endif
            //for (iDim=0;iDim<NDIM;++iDim)
            //for (iCDim=0;iCDim<NDIM;++iCDim){
            //    if ((node+iLNode)->ja[iDim][iCDim]<TINYNUM)(node+iLNode)->ja[iDim][iCDim]=0;
            //    if ((node+iLNode)->J<TINYNUM) (node+iLNode)->J=0;
            //}
            iType=0;
            for (iDim=0;iDim<NDIM;++iDim){
                dl[iDim]=0;
                iType+=abs((node+iLNode)->locTyp[iDim]-BOUNSTE);
            }
            for (iDim=0;iDim<NDIM;++iDim){
                //if (iType==0){
                    for (i=0;i<NDIM;++i) dl[iDim]+=invJ[i][iDim]*invJ[i][iDim];
                    dl[iDim]=sqrt(dl[iDim]);
                    dlMax=max(dlMax,dl[iDim]);
                    dlMin=min(dlMin,dl[iDim]);
                //}
            }
        }
    }
    #ifdef _MPI
    MPI_Allreduce(&dlMax,&dlGlMax,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);
    MPI_Allreduce(&dlMin,&dlGlMin,1,MPI_DOUBLE,MPI_MIN,MPI_COMM_WORLD);
    #else
    dlGlMax=dlMax;
    dlGlMin=dlMin;
    #endif
    m_msg="Maximum mesh size is "+to_string(dlGlMax)+" and minimum mesh size is "+to_string(dlGlMin);
    logCtrl.info(m_msg.c_str());
    if (inputCtrl.dt>0){
        m_msg="dt = "+to_string(inputCtrl.dt)+"  is given in control file";
        logCtrl.info(m_msg.c_str());
        inputCtrl.CFL=(inputCtrl.Ma_ic+1)*inputCtrl.dt/dlGlMin;
        m_msg="CFL = "+to_string(inputCtrl.CFL)+"  is calculated by dt";
        logCtrl.info(m_msg.c_str());  
    }else{
        inputCtrl.dt=dlGlMin*inputCtrl.CFL/(1+inputCtrl.Ma_ic);
        m_msg="dt = "+to_string(inputCtrl.dt)+"  is calculated based on CFL ";
        logCtrl.info(m_msg.c_str());
    }
    m_msg="Peroid time = "+to_string(1.0/inputCtrl.srcFreq)+" Steps = "+to_string(1.0/inputCtrl.srcFreq/inputCtrl.dt);
    logCtrl.info(m_msg.c_str());
}